Mott insulator dynamics 
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The hydrodynamics of a lattice Bose gas in a time-dependent external potential is studied in 
a mean-field approximation. The conditions under which a Mott insulating region can melt, and 
the local density adjust to the new potential, are determined. In the case of a suddenly switched 
potential, it is found that the Mott insulator stays insulating and the density will not adjust if the 
switch is too abrupt. This comes about because too rapid currents result in Bloch oscillation- type 
current reversals. For a stirrer moved through a Mott insulating cloud, it is seen that only if the 
stirrer starts in a superfiuid region and the velocity is comparable to the time scale set by the 
tunneling, will the Mott insulator be affected. 



I. INTRODUCTION 

When Greiner et al. put a system of bosonic atoms 
in an optical lattice and demonstrated the Mott transi- 
tion research on the physics of strongly correlated 
systems entered a new phase. Cold atoms in optical lat- 
tices offer control and detection techniques unthinkable 
in a condensed matter system. The versatility of optical 
lattices have teased the imagination of theorists, and co- 
pious amounts of papers suggesting exotic types of quan- 
tum phase transitions have been published in the past ten 
years jH. However, almost all of the suggested quantum 
phases have proven impossible to realize experimentally 
at the present for various technical reasons. For example, 
many rely on spin ordering, which requires temperatures 
far below what is possible at the present (see, however, 
the ingenious experiment by Simon et al. using occupa- 
tion as an effective spin degree of freedom to realize a 
quantum Ising model [1]). 

While the quest for achieving lower temperatures is on 
Q , there is still lots of interesting science to be done with 
the relatively simple Bose-Hubbard model on a square 
lattice. With the possibilities of monitoring directly the 
local density and low-order correlation functions on the 
one hand [3, Q , and the possibility to manipulate exter- 
nal potentials in real time on the other hand , the idea 
to study the real-time dynamics of a strongly correlated 
system suggests itself. In particular, these systems are by 
default trapped in inhomogeneous potentials and there- 
fore they contain finite regions of superfiuid and Mott in- 
sulator sitting side by side. With current techniques, one 
can investigate how the interfaces between such regions 
behave in real time under various kinds of perturbation. 

The majority of studies of dynamics of la,ttice Bose sys- 
tems have been concerned with quenches [8|, Bloch oscil- 
lations , and collective modes in traps [10l - [l2| . These 
aspects of trapped lattice bosons are by now well un- 
derstood. Worth mentioning is also an important study 
of the critical current in interacting lattice Bose systems 
[isj . Concerning the macroscopic transport of matter 
and redistribution between quantum phases due to po- 
tentials and currents - i.e., the hydrodynamics - the lit- 
erature is less complete. Natu et al. study the redistribu- 



tion of matter within an optical lattice following a quench 
. Fischer et al. were concerned with the velocity of a 
moving superfiuid- Mott insulator interface [15| . Karlsson 
et al. study the behavior of a trapped partly Mott insu- 
lating system after turn-off of the trap 16], and Snock 
and Hofstetter study the dynamics upon displacement of 
the trap this paper will be further discussed below. 

The literature survey above indicates that there is need 
for a more exhaustive understanding of how Mott insu- 
lators react to macroscopic currents and potential gradi- 
ents, which is the subject of this paper. It will provide 
two sets of numerical examples of how a Bose-Hubbard 
system containing coexisting Mott and superfiuid regions 
evolves in time following a perturbation in the potential. 
In order to do this, we evolve the Bose-Hubbard Hamil- 
tonian in time within the mean-field Gutzwiller approxi- 
mation. In Sec.|TT]the governing equations of motion are 
put up. In Sec. IIIIl I study the in- and outflow to or 
from a Mott insulator following the onset of a potential 
gradient. Sec. |IV] discusses the perturbation of a Mott 
insulator using a localized stirrer. In Sec.|V]l summarize 
and conclude. 



II. LATTICE BOSE GAS 

We study a system of zero-temperature bosons hopping 
on a lattice and in addition subject to a more slowly vary- 
ing external potential. The many-body Bose-Hubbard 
Hamiltonian is 
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Here, J is the tunneling matrix element, U is the inter- 
action parameter and /i is the chemical potential. V{rj) 
is the external potential and the index j runs over the 
lattice sites, rj is the spatial coordinate of the j:th lat- 
tice site. The sum subscripted < > runs over nearest 
neighbors. We work in units in which h — 1 and the lat- 
tice constant is also unity; hence, both frequencies and 
velocities can be measured in units of U. This paper 
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will be concerned with a two-dimensional square lattice. 
However, note that the mean- field approximation that 
will be used here is most accurate in higher dimensions, 
so the results obtained are best seen as qualitatively de- 
scribing 3D physics. 

The Gutzwiller approximation is based on a mean- field 
ansatz for the many-body state 

i^GW)=ni'/'.w)- (2) 

3 

In the numerical computations, the on-site states are ex- 
panded in a local Fock basis with an upper cutoff nmax, 

■"max 

l-^.W) - E c,At)W)r- (3) 

In order to calculate the ground state, the Hamiltonian 
([1]) is minimized with respect to the complex coefficients 
Cj^n- The time development is obtained by pro pagating 
the coupled equations of motion [l^, [13, [l9l - [2l| 

BC r) 

^^-gC^i^Gimii^Git)). (4) 

In order to diagnose the state, the local total density rij 
and condensate wave function V&j = (dj) are computed. 
We will characterize the system by studying the behavior 
ofrij, the local condensate density ricj = l^'j and phase 
if — arg^. Although simplistic and certainly quantita- 
tively inaccurate in low dimensions, I see this mean- field 
study of the time development as a first study which can 
later be vindicated or falsified in more accurate simula- 
tion schemes. 

The phase diagram for the Bose-Hubbard model was 
discussed in, e.g., Ref. (2^. The mean- field critical point 
at n = 1 in 2D is t « 0.042f/ at ^ = 0.5U. If the sys- 
tem is in an external potential, the density is usually well 
modeled by a local-density approximation (LDA), so that 
alternating superfluid and Mott insulating regions exist 
alongside each other, determined by the local chemical 
potential fi — V{r). The density profile assumes a char- 
acteristic wedding-cake structure [23| . 

III. OUTFLOW FROM A MOTT INSULATOR 

We consider a 2D system trapped in a harmonic po- 
tential, 

^o(r,) = y(:.|+yf), (5) 

where the integer coordinates Xj and i/j run from —L + 1 
to L as j runs from to - 1. With the choice J = 
0.03J7 and n = 0.5U, the system is Mott insulating in the 
center with a surrounding superfluid shell. The initial 
condensate density profile, with L = 32 and uj = 0.07J7, 
is graphed in Figure [IJa). At time t — 0, with the system 



(a) (c) 




-20 20 



X 



(b) (d) 




200 400 600 800 
t 



FIG. 1. [Color online] Time development of a 2D lattice Bose 
gas following the turn-on of a wide Gaussian perturbation po- 
tential, (a) Spatial distribution of the condensate density be- 
fore the switch of the potential, (b) Final condensate density, 
after an evolution of duration t — 800U~^ . (c) Cross-section 
of the final total density profile (symbols) and potential pro- 
file after the switch (line), (b) Fraction of Bose-Einstein con- 
densed atoms as a function of time. The tunneling matrix 
element is J = 0.03!7, the trap frequency is a; = 0.07(7, the 
chemical potential is /x = 0.5J7, and the Gaussian potential 
has width W = 10 lattice sites and strength Vi = 0.5(7. 

in the ground state of the trap, an additional perturbing 
potential is suddenly switched on, 

V{r,) = Voirj) + yie-("?+^?)/^'. (6) 

The resulting total potential assumes a toroidal form, as 
can be seen in Fig. [Ijc). In the case of Fig. [T] I chose 
Vi = 0.5(7 and W = 10, to make a relatively wide per- 
turbation. There results a rather violent time develop- 
ment, but over time it is seen that the Mott insulating 
phase is molten and the system assumes something rem- 
iniscent of a steady state; Figs. [Ub)-(c) indicate that 
apart from the noise induced (which I would like to call 
thermal noise, although the relation between the current 
mean-field approximation and a finite-temperature Hub- 
bard model is not at all clear) the final density follows 
the potential profile, closely approximating the LDA. In 
Fig. [TJ^d) , the total condensate fraction Nc = ricj is 
plotted as a function of time. As expected it is increas- 
ing as the Mott insulator melts, and the time scale of the 
process is of the order of 100?7~^ or 10J~^. 

However, the picture is different if the potential is 
switched to large enough values. In Fig. [2l the final po- 
tential strength is now set to the larger value Vi — 1.0(7. 
In this case, the Mott insulator resists melting and the 
system is hindered from approaching equilibrium. A 
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FIG. 2. [Color online] Time development of a 2D lattice Bose 
gas following the turn-on of a wide Gaussian perturbation 
potential. Panels and parameters are as in Fig. [U except the 
height of the Gaussian potential, Vi — l.OU. 



steady state seems to be attained, again after a few in- 
verted hopping periods J~^, but it is not the LDA distri- 
bution seen in the previous numerical experiment. Figure 
[2ljb-c) show that there is still left a Mott insulating re- 
gion with unit density at the trap center. In Figure HJd) 
we see that in this case, the condensate fraction has also 
increased a bit, but the important finding is that the cen- 
tral Mott insulator is not entirely molten. This behavior 
persists over a range of values of the inverted final po- 
tential. Further numerical experimentation shows that 
when Vi is larger than about 0.8U, the Mott insulator 
insulates. 

We have thus found that if the potential is changed a 
little, the Mott insulator will melt, but if it is changed 
a lot, it will not. The reason for this behavior becomes 
clearer if one studies the velocity of the gas. In Fig. [3] 
the phase of the condensate is plotted for the first nu- 
merical experiment. In Fig. |4] it is plotted for the sec- 
ond one. The fluid velocity is related to the wavenum- 
ber of the phase variation pattern through the relation 
Vj — Jsinkj, where j = x,y denotes a Cartesian coordi- 
nate axis. The wavenumbers kj are restricted to lie in the 
interval (— tt, tt). If the potential is steep, the fluid will be 
accelerated to the maximum wavenumber tt, as is seen to 
be the case in Fig. |4l This gives rise to a current reversal, 
which is in essence a Bloch oscillation. At longer times 
turbulent processes arise, possibly aided by the dynam- 
ical instability known to take place in a condensate at 
k — Ti/2 0], and in interacting systems at even smaller 
wavevectors p^ . The subsequent time development is 
noisy. This in itself would not necessarily stop the fluid 
from eventually flowing into the newly created potential 
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FIG. 3. [Color online] Spatial distribution of the condensate 
phase of a 2D lattice Bose gas in a harmonic trap plus Gaus- 
sian potential. Parameters are as in Fig. [U featuring the 
weaker Gaussian potential with Vi = 0.5!7, and the panels 
refer to times (a) t = 20I7"\ (b) t = 40)7"\ (c) t = 60(7"\ 
and (d) t = lOOC/"^ 
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FIG. 4. [Color online] Spatial distribution of the condensate 
phase of a 2D lattice Bose gas in a harmonic trap plus Gaus- 
sian potential. Parameters are as in Fig. [21 featuring the 
stronger Gaussian potential with V\ = 1.0(7, and the panels 
refer to times (a) t = 20(7"\ (b) t = 40(7"\ (c) t = 60f/"\ 
and (d) t = 100f/-\ 



wells, and it does not do so if the system is purely super- 
fluid (as seen in simulations not shown here). However, 
when a Mott insulating region is in the way, the simu- 
lations indicate that the whole process is halted and the 
Mott insulator insulates. 

Snoek and Hofstetter [17| studied the dynamics of a 
trapped lattice boson system following lateral displace- 
ment of the trap. Similarly to the present study, they 
found that Bloch oscillations play a role for the dynam- 
ics in the case of large displacements. However, in that 
study, the system always relaxed towards the new equi- 
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FIG. 5. Final central condensate density after the turn-on of 
a wide Gaussian perturbation potential and subsequent evo- 
lution for a duration of 800C/~^. The condensate density Uc 
is averaged over the nine centermost points of the lattice. 
Plusses denote the result of simulations with a Gaussian per- 
turbation with width W = 10 and height Vi; asterisks are for 
a Gaussian with width W = 15, and circles are for a Loren- 
zian profile with width W = 15. In (a), the x axis is the 
height of the perturbing potential, Vi; in (b), it denotes the 
maximum potential slope, V^^x- 
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FIG. 6. [Color online] Time development of a 2D lattice Bose 
gas following the turn-on of a wide Gaussian perturbation 
potential. Panels and parameters are as in Fig. [1] except J = 
0.02(7, fi = 1.5U, and the height of the Gaussian potential is 
Vi = 1.0(7. 



librium in two dimensions. In the present setting, we see 
that this behavior can be violated. 

A wider scan over parameters is collected in Fig. [S] 
Here, we have measured the condensate density Uc in the 
center at the end of the simulation, and plotted it against 
the potential height Vq as well as against the maximum 
slope T^ax- Two different Gaussian perturbations are 
chosen, with widths W — 10 and W — 15. In addition, 
I tried a perturbation with a Lorenzian profile, not be- 
cause of experimental relevance (it is probably very hard 
to make), but in order to show that the qualitative re- 
sult is insensitive to the shape of the potential. It is 
seen that the core first becomes superfluid when the po- 
tential height Vi exceeds about OAU, simply because of 
LDA considerations. Then, for large enough values of Vi, 
there is no superfluid in the core, as we have seen above. 
This cutoff value depends on the specific potential, as 
seen in Fig. [SJa), but in Fig. [Sib), the curves collapse 
onto each other when the condensate density is instead 
plotted against the maximum potential slope, V^max: i-^-j 
the maximum potential difference between adjacent sites. 
The critical value of Vj^^x this case approximately 
0.08J7, not an integer multiple of U, J, nor u, ruling out 
the simplest guesses for resonance physics, but consistent 
with the picture that accelerated atoms experience Bloch 
oscillations. 

So far we have only studied a single combination of 
tunneling strength J, trap frequency uj and chemical po- 
tential /i, giving a Mott insulator with n = 1. Starting in 
a different Mott insulating region gives the same result. 



as seen in Fig. |6l Here, the parameters are chosen as 
J = 0.02C7, = 1.5U, and Vi = l.OU. This means that 
the central region in the initial state is a n = 2 Mott insu- 
lator surrounded by three shells of alternating superfluid 
and Mott insulator. Numerical experimentation shows 
the same behavior as reported above: For a weak enough 
potential, the central Mott insulator is molten, but for a 
stronger one it is not. 

More curious behavior is seen in a system with a super- 
fluid region in the center. We choose J — 0.3U, /i = l.OU, 
and 14 = 1.5U to produce Fig. [T] As the peaked potential 
rises in the center, bosons are flowing out from the cen- 
tral superfluid and flnally into the outermost superfluid 
shell. However, a partly Mott insulating region is left in 
the center, with a few superfluid atoms intermixed, and 
then, subsequent dynamics is halted. It is seen in Fig. [7] 
that the final steady state does not rhyme very well with 
the applied potential. A closer look at the time depen- 
dence is given in Fig. [51 It is seen that the superfluid 
atoms flow out by creating four jets through the Mott 
insulator, following the spatial symmetry of the lattice. 
A noisy state is left behind, containing a small amount 
of condensate, but as seen above, the final steady state 
clearly violates the LDA. In this case, it appears that 
the outflowing superfluid bosons leave a Mott insulator 
behind; however, a quantitative explanation seems to be 
out of reach at the present. 
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FIG. 7. [Color online] Time development of a 2D lattice Bose 
gas following the turn-on of a wide Gaussian perturbation 
potential. Panels and parameters are as in Fig. [T] except J = 
0.02(7, fi = l.OU, and the height of the Gaussian potential is 
Vi = 1.5U. 




FIG. 9. [Color online] Time development of a 2D lattice Bose 
gas as a Gaussian "spoon" potential is moved through it. The 
potential has width W = 3, amplitude Vi — 5U, and velocity 
V = 0.517. The remaining parameters are as in Fig. [T] 




FIG. 8. [Color online] Spatial distribution of the condensate 
density nc of a 2D lattice Bose gas in a harmonic trap plus 
Gaussian potential. Parameters are as in Fig. [71 featuring 
several superfluid and Mott insulating shells, and the panels 
refer to times (a) t = 0(7" \ (b) t = 30(7" \ (c) t = 60(7" \ 
(d) t = 90(7"\ (e) t = 120(7"\ and (f) t = 600(7"^ 



IV. STIRRING A MOTT INSULATOR 

Now let us investigate what it takes to make a hole in 
a Mott insulator. To this end, we take a system with fi = 
0.5 as above; as we have seen this gives a Mott insulator 
in the center with an approximate width of 20 lattice 
sites. We move a narrow Gaussian "spoon" potential 
through the cloud, as was done in Ref. [2^ to excite 
vortex pairs in a BEG (cf. To produce a reasonably 

narrow and strong stirrer I choose the parameters W = 3 
and Vi = 5.QU, and the center of the potential moves 
linearly through the cloud with a velocity v, 

V{rj,t) = Vo{rj) + yie-«"^-''*)'+^?)/^'. (7) 

First, we conclude that a Mott insulator will not be per- 
turbed unless the spoon starts in a superfluid region. 
This is in accordance with the definition of Mott insu- 
lator, and a series of numerical experiments (not shown 
here) have confirmed this. Thus, we start the spoon at 
the edge of the simulation cell and move it through the 
superfiuid shell and then to the Mott insulating interior 
of the cloud. 

Figure [9] shows what happens when the spoon is moved 
quickly through the cloud: The Mott insulator insulates, 
and the presence of the spoon is only felt as it passes 
through the thin superfiuid shell. In this simulation we 
chose a velocity v — 0.5U (recall that h = 1 and the 
lattice constant is 1). This is, in fact, too fast for the 
tunneling dynamics to respond; in order to melt the in- 
sulator one needs to move across a lattice site at a time 
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FIG. 10. [Color online] As in Fig. (9] but with velocity v = 
0.025[/. 
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FIG. 11. Fraction of Bose-Einstein condensed atoms at the 
end of a sweep with a Gaussian "spoon" potential, as a func- 
tion of the velocity with which the spoon is moved. Dashed 
line: Initial condensate fraction. 

scale comparable with J^^. Such a case is shown in Fig- 
ure [TUl where v = 0.025U. Here we see how the spoon, 
given enough time, excites bosons out of the Mott insula- 
tor and depletes the density in the vicinity. An attempt 
to quantify this result is made in Fig. 1111 Here, we mea- 
sure the fraction of condensed atoms Nc/N at the final 
time, when the spoon has traversed the entire simulation 
cell with the length of 64 sites. It is seen that the upper 
critical velocity for exciting atoms into the condensate is 



around v = O.IU, which is comparable to the tunneling 
J. We note that the sound velocity is likewise of order J. 
These simulation results indicate that making a hole in 
the Mott insulator is basically a question of moving slow 
enough, or else the Mott insulator will not budge. When 
the velocity is decreased below the tunneling strength, 
though, the curve is seen to decrease again. Indeed, it is 
natural to expect that for a slow enough spoon, the fluid 
will adjust adiabatically and the perturbation will again 
be minimized. Thus, in the case of a moving spoon, we 
have identified three parameter regimes: An adiabatic 
regime for velocities about an order of magnitude below 
J; a supersonic regime in which the edge superfluid has 
no time to react; and an intermediate regime in which 
the spoon is slow enough to excite the Mott insulating 
atoms and fast enough to do maximum damage. 



V. CONCLUSION 

In summary, I have studied the hydrodynamic behavior 
of a system of lattice bosons in an external potential, by 
means of a number of numerical experiments. In the first 
set of experiments, the time development is monitored 
after a sudden switch of the potential. It is found that if 
the potential is switched by a small enough amount, the 
bosons will move in order to adjust to the new potential 
and a Mott insulator may melt. However, if the switch is 
large enough, the Mott insulator will not give in and the 
system stays in a metastable configuration. The reason 
for this behavior is that the supercurrent breaks down 
in a Bloch oscillation if the potential switch is too large. 
In the second set of numerical experiments, a "spoon" - 
type of localized potential is moved through a system 
containing a Mott insulating region. The Mott insulator 
is affected by the spoon if is moved at just the right 
pace, so that the time scale for traversing a lattice site is 
comparable to the tunneling time scale, J~^. For faster 
spoons, its presence is not felt in the Mott insulator, and 
for slower spoons, the whole process is adiabatic. 

Experimentally, present techniques for in-situ measure- 
ment of filling factor Q is the most straightforward 
way of testing the predictions made here. A classic time- 
of-flight measurement, which in effect can tell the Bose- 
Einstein condensed fraction of atoms Ij, is also feasible, 
as can be seen from the plots of said quantity in the fig- 
ures. 
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